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I. INTRODUCTION 

The fully quantized Rabi model [1, 2] constitutes probably the simplest strongly coupled quantum system with an 
infinite dimensional Hilbert space. It is described by the Hamiltonian (h = 1), 

H R = uia^a + ga x (a + a t ) + Aa z . (I) 

Hr possesses a Z 2 -symmetry which renders it integrable [3] , although the associated two invariant subspaces (parity- 
chains) are each infinite-dimensional themselves. The regular spectrum of (1) is determined in each parity chain by 
the zeroes of a transcendental function G± (x) , which can be expressed through confluent Heun functions [4] . G± (x) 
has poles at integer multiples of u, which allows to obtain results on the distribution of eigenvalues and the location 
of degeneracies which are related to the quasi-exact (exceptional) spectrum [3] . The approach based on the symmetry 
of the model leads in this way to a unified picture of all spectral properties. 

Independent from the qualitative understanding of the physics governing the quantum Rabi model is the question 
whether the analytical solution could invalidate the widely used numerical determination of the spectrum through 
exact diagonalization on a truncated state space, or the equivalent continued fraction techniques [5-12]. There has 
been recent confusion on this point [13, 14] and the present paper intends to clarify the situation regarding the 
methods based on continued fractions. These fall into two classes, which we denote with A and B. 

II. METHOD A 

One considers the full model, i.e. the state space is assumed to be infinite-dimensional, W = L 2 (R) <E) C 2 . Work 
in this direction goes back to the seminal paper by Schweber [5]. He used the isomorphism between L 2 (R) and 
the Bargmann space B of analytic functions [15] to obtain a differential equation for the spin-up component of the 
eigenfunction with energy E, 

oo 

j, E ( z ) = e -9*/» K n {E)(z + ^) n , (2) 

n— — oo 

leading to a linear three-term recurrence relation for the coefficients K n (E), 

nK n = f n -!(E)K n -i - K n - 2 , (3) 

with 

f n (E) = ^ + ±(nuj-x(E) + —^ ) (4) 

uj 2g \ x(E) — nuj J 

and x(E) = E + g 2 /uj. The same equation results if the Z 2 -invariance is used to project the Hamiltonian to one of its 
invariant subspaces with fixed parity [3]. Each of the parity-chains T-L± is isomorphic to B [3] and we set H = H + = B 
in the following, confining the discussion to positive parity. The energy eigenvalues E n , n = 0, 1,2, . . . of Hr are 
determined by imposing on iPe(z) the condition to be an element of B, henceforth called the "spectral condition". 
This condition consists of two parts, the first is the familiar normalizability, |V>.e(z)| < oo, where | • | denotes the norm 
in B. Because the differential equation for "4>e(z) in (2) has an irregular singular point of rank one at infinity [16], 
4>e( z ) ~ e cz z p (c + Ci/z + . . .) with complex constants c, p, c , C\ [17]. As B contains all functions growing like e cz for 
\z\ — > oo, iPe(z) fulfills the first part of the spectral condition for any real E. The second part demands iPe(z) to be 
analytic everywhere in C, and it is this second part which fixes the discrete set {E n } ne iN Q , the spectrum of Hr. The 
differential equation for i/je(z) has two regular singular points at z = ±g/oj and iI>e{z) will be analytic in C if it is 
analytic at both of these points. iPe(z) is expanded in (2) around z = —g/uj and analytic there if K n = for n < — 1. 
This leads to the initial condition for the recurrence (3): Kq = 1 and K\{E) = fo(E). It yields still no equation for 
E, which is obtained by demanding analyticity at the second singular point z — g/u>. From the Poincare analysis 
of (3) it follows that there are two possibilities for the limiting behavior of the K n (E): K n +\(E) / K n (E) — > u)/(2g) 
or K n+ i(E)/K n (E) — > 0. The latter behavior corresponds to the unique minimal solution {K™ ln } ne ^ of (3) with 
initial condition Ko — 1 [18]. Note that K™ m is determined in terms of {fi(E), /2(E), . . .} for any E and does not 
depend on fo(E) [19, 20]. All other solutions {K n } n ^ of (3) are dominant, i.e. Xm\ n ^ toa K™ va jK n = 0. In general, 
K™ m (E) f (E), but for E G {E , E x , . . .} one has K™ m (E n ) = f (E n ). This condition determines the spectrum 
because only then ipE n (z) is analytic at both —g/uj and g/u and therefore in all of C. Analyticity at —g/u> fixes 
K\ = fo(E) and analyticity at g/ui enforces an infinite radius of convergence of the expansion (2), which would be 
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finite (R — 2g/uj) for any dominant solution with K n+ i(E) / K n (E) — > u)/(2g). Only the minimal solution yields 
R = oo, therefore Ki(E) = K™ m (E) = fo(E), which is the sought after equation for E. 

To compute the minimal solution of (3), Schweber proceeds to represent the quotient K\(E)/Kq through the infinite 
continued fraction 

= f (E) = l\h(E) - 2\f 2 (E) - 3\f 3 (E) - . . . . (5) 

As pointed out in [3], this equation amounts to a tautology if the formal expression on its right hand side is not 
augmented with a prescription for the asymptotics of the tail £,n(E) in 

Foo(E) = l\h(E) - 2\f 2 (E) - 3\f 3 {E) - . . . = l\h(E) - 2\f 2 (E) - ... (N - l)\f N -i(E) - £ N {E)- (6) 

Fortunately, the standard definition of an infinite continued fraction (E) as the limit (if it exists) of a sequence of 
finite continued fractions F n (E), 

Foc(E) = lim (l\.h(E) - 2\f 2 (E) n\f n (E)) = lim F n (E) (7) 

yields the correct prescription. The reason is that the limit (7) exists in case of the asymptotic behavior 

^T^S=° (8) 

but because = nK n /K n _i and f n ~ nui/(2g) for large n, this can only happen if lim^oo K n {E) / K n _\(E) = 0, 
which characterizes the minimal solution. Pincherle's theorem, which connects minimal solutions with Kq ^ to 
continued fractions is a simple consequence of the uniqueness of {K™ m } [19]. In other words, the continued fraction 
(7) is a short-hand notation for the minimal solution K ™ m (E)/Kq of the recurrence (3), which is defined for arbitrary 
E. Schweber has proven the (pointwise) convergence of (7) using a corollary of Pringsheim's main theorem [21, 22]. 
However, in contrast to power series expansions, the actual value of (7) can only be obtained by checking convergence 
numerically for each E separately, by computing a set of F n (E) with large enough n. It is not possible to deduce 
from the smallness of the tail 



UE) = n\f n (E) - (n + - • . • (9) 

the value of Foo(E) for given F n -i(E) to a prescribed accuracy because the continued fraction itself may diverge 
while stays finite [23]. Therefore, the actual computation concerns always some F N (E) with finite N. It is easy to 
see that the equation fo(E) = Fn(E) is just the spectral condition for the Hamiltonian Hn = PnFIrPn, where Pn 
denotes the projector onto the finite-dimensional Hilbert space Hn, spanned by functions of the form 

^(z) = e- BZ ^(z+^Y , n = 0,1,2, ... N. (10) 

If ?Pe(z) = X^^Lo Kn4>n{z) satisfies formally (Hr — E)iPe{z) = 0, it follows 

oo oo 

K n H R( j> n = [(2.9/n(£) + E)K n - 2g[(n + l)K n+1 + K n ^}] <j> n , (11) 

n=0 n=0 

because of (3). Obviously, PnHrPn^n+i — 0, which entails that the coefficient of <pN in (11) becomes (2gfx(E) + 
E)Kn — 2gKjy_i, which means K^+i = 0, i.e. the recurrence (3) is cut off at Kn- Now (3) together with the initial 
condition K\ = fo(E) and additionally Kn+i = is equivalent to an equation involving a finite continued fraction, 

f (E) = l|/x(£7) - 2\f 2 (E) - 3\f 3 (E) - . ..N\f N (E) = F N (E), (12) 

and provides therefore the spectral condition in the space Hn- 

Note that this space does not correspond to the span of a finite subset of Fock states , n — . . . N} as the 
truncated spaces HS NS> used in method B, but all elements of Hn can be approximated to arbitrary precision by 
elements of VS N ^ for sufficiently large N (see below). Each solution E^ of (12) approximates some solution E n of 
F rx (E) = fo(E) to precision e if TV > N{n,e). The latter property follows from the fact that the convergence of 
the tail £n(E) depends on the value of E being smaller than a certain upper bound depending on N (see [5] or 
Eq. (24) below). The solution E^ of (12) must be smaller than this bound for Fn(E^) to be a good approximation 
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of F OQ (E^). However, for each finite E n being a solution of the exact equation fo(E) = F OQ (E), there is a N(E n ) 
such that for all N > N(E n ) the function F N (E n ) approximates F 00 (E n ) to a given precision and E n rts E^ will be 
thus close to the solution E% of (12). 

One may formalize this statement as follows. The domain of the operators Hn is extended to % by setting 
ELn^tu = for m > N + 1. Then the spectrum (j{Hr) = {E n } ne ^ of the operator Hr is the limit of the spectra 
<t(Hn) = {En}neiN in the sense that there exist a number M(N), unbounded in N, such that 

lim max \E* - E n \ = 0. (13) 

JV-K» n<M(JV) 



The spectral condition on "Hat provides an approximation to the spectrum of the full model precisely because it 
approximates the exact condition of analyticity in the Bargmann space. It follows that the continued fraction approach 
initiated by Schweber and pursued by several authors [7-11] can be used to obtain the spectrum of the quantum Rabi 
model to arbitrary precision. The method does not use the Z2-symmetry of Hr as the Bargmann condition is 
implemented via (5). Indeed, (5) only depends on A 2 and does not discern the two parity chains, which is the reason 
why the relation of isolated (Juddian) solutions at x(E) = nui to level crossings, i.e. reducible representations of Z2, 
appears to be accidental within Schwebers approach [24] . The infinite continued fraction Foo (E) is completely opaque 
to an analysis of its general behavior as function of E: F QO (E) is the formal quotient of two functions, 

Foo (E) = lim £g| (14) 

n-s-oo t> n {E) 

where A n , B n both satisfy the recurrence relation 



C n = /„(S)C n _i - nC„_ 2 , n > 2 (15) 

and j4 = 0,^4_i = l,B = 1,B_ 1 = 0. Only the quotient of A n and B n is well-defined for n — > 00. Therefore the 
location of zeros and poles of F oc (E) — fo(E) cannot be read off from the behavior of numerator or denominator 
in the limit n — > 00, in contrast to the function G±(x) whose pole structure is known and leads to simple rules for 
the distribution of eigenvalues. Method A can be used to calculate the spectrum of the Rabi model to any desired 
accuracy, but provides no insight beyond direct numerically exact diagonalization in a finite-dimensional state space 
with dimension N + 1, if (5) is truncated at the order N. 



III. METHOD B 



This approach goes back to Swain [6]. It is based on a three-term recurrence relation for the determinants of 
tridiagonal matrices. Because the quotient of two such determinants can be interpreted as a matrix element of the 
resolvent of H R N \ if defined on a truncated Hilbcrt space, H\ N ^ = C N+1 ® C 2 , one may write this matrix element as 
a finite continued fraction. By repeating Swain's calculation, the authors of recent work [12, 14] tie the "solvability" 
of the Rabi model to his representation of the resolvent. It is further claimed in [14] that this representation, which is 
defined in %[ N ^ only, has a well-defined limit N — > 00 at least for g 2 < lo and provides in this way an exact solution 
of the Rabi model without making use of the analyticity condition employed by method A. This claim is based on the 
belief that the operator Hr, defined on an infinite-dimensional Hilbert space, can be approximated by the sequence of 
finite-dimensional matrices H R , which is incorrect because Hr is not compact [25]. Nevertheless, the full spectrum 

of Hr can be obtained by method B as well — but this is entirely due to the equivalence of the spaces %\ to the 
spaces %n <£> C 2 from method A — and rests therefore again on the Bargmann criterion, as shown in the following. 

There are several equivalent formulations of Swain's approach [6, 12, 14, 26]. In its simplest form [12], one considers 
the Hamiltonian in each parity chain separately, where it corresponds to a tridiagonal matrix. Define the parity chain 
%± = span{|n, ±(— l)™)}„ e iN with \n, 1) = \n) ® |t) and \n, — 1) = \n) ® \\). The Fock states \n) are eigenstates 
of a^a: a) a\n) = n\n). %± is isomorphic to i 2 (K), i.e. to the span of Fock states \n) for n € IN. To define 

the finite-dimensional approximants to Hr in %±, one uses projection operators P±*\ projects onto 

U { ± ] - span{|n,±(-l)")} < „<jv- The truncated parity chain is then isomorphic to C N+1 and the truncated 
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Hamiltonian in H^ reads H { ± ] = Pi N) H R Pi N) . 



assumes the tridiagonal form 



/ bf v«T 



Hi N) = M± = 



i_ 

C12 



N — l 



\ 



/OLN 



'N 



(16) 



with 



&f =jw±(-l)'A, a^jg 2 , j =0,1, 



. AT. 



(17) 



To derive a recurrence relation for det Mq , define matrices by deleting the first j rows and columns from Mq . 
Then it follows, 



dctM; 



bf dct Mf +1 - a J+1 det Mf +2 



(18) 



for j = 0, . 



.N, setting detM± +1 - ^,, UAr+2 



1, detM^ ,„ — 0. This three-term recurrence can be turned into a non-linear 



two-term recurrence for G- 



AetMf +1 /detMf, 



G, 



bf - a j+1 Gf +1 



(19) 



To compute the resolvent of H { ±\ we set Mq' = E - H^> and find for Gf(E) = (0,±1\(E - H^'y^O, ±1) a 
representation in terms of a finite continued fraction 

G±(E) = l\b±(E) - ai \bf(E) - a 2 \bf(E) a N \b%(E), (20) 



(JV) 



where the b • (S) , aj are now defined as 



bf(E)=E-jwT(-l) j A, 



39 



j = 0,l,...N. 



(21) 



We drop the parity index from now on and write Gq N \e) for the expression (20). Gq N \e) = det M\j detM has 
poles at the eigenvalues of H^ N \ provided that the pole (a zero of det(£ - ffW) at E n ) is not lifted by a zero of 
det Mi, which could occur if (0, ±l\ip n ) — for the eigenstate \ip n ). This does not happen in the present case because 
all dj are nonzero. Nevertheless, the numerical evaluation of (20) is compromised by the rapid decay of the overlap 
(0,±l|^ n ) with growing n. It becomes apparent in Fig. 1 of [14], where the peaks of the spectral density are very 
small already for n = 9, 10. 

Method B presumes that the poles of limjv^oo Gq N ^ give the spectrum of the untruncated Rabi model if G^(E) = 
lini/v-j.oo can be shown to exist for any real E. To prove this, Ziegler [14] invokes the main theorem of Pringsheim 
[22] , which states that the sequence of finite continued fractions 



F N = ai\fii - a 2 \/3 2 - a 3 \j3 3 



■ OtN\PN 



(22) 



converges for N — > oo to a value F with \F\ < 1 if \fij\ > \otj\ + 1 for all j = 1,2, 
dimcnsionless, this theorem could be used to compute the tail £ r 



If the a,' 



of Gg° 



Pj would be 



*n+l|"n+l 



for sufficiently large n and fixed E, because \bj \ ~ jco and \a,j 



jg 2 for j > n » 1 such that \bj\ 



(23) 

> laA + 1 for all 



j > n, provided g 2 < oj. The convergence of £ n {E) then ensures convergence of the full continued fraction G^{E) in 
the generalized sense, which is sufficient to conclude that the distribution of the poles of G^(E) is well-defined. 

To obtain dimensionless quantities, g and A in (1) are scaled with co: g — g/uj, A = A/w, u> = 1. The condition 
g 2 < oj turns into g 2 < 1, or \g\ < uj. This condition seems to invalidate the proof for the deep strong coupling regime 
\g\ > ui, which is of current interest [27]. But Ziegler [14] does not rescale the Hamiltonian and obtains instead of 
\g\ < uj the meaningless expression g 2 < ui, which compares quantities of different dimension. However, it is easy to 
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correct this mistake by using the following generalization of Pringsheim's theorem for dimensionful quantities which 
lifts at the same time the restriction \g\ < lu [28]. 

Theorem: Consider the continued fractions (22) and assume for the dimension [Fn] of Fn the relation [Fn] = \(3j] 
and [a-,] = [F N ]. If IIJ=i \ a j\ 1S unbounded in n and the inequalities |/3j| > |aj|/c + c hold for all j > 1 and some 
constant c with dimension [c] = [/3j], then liniAr^oo Fn = F exists and |F| < c. 

The theorem is applicable to the tail £ n (F) of Gq'(E). Moreover, the possibility to choose c freely leads to the 
following lower bound for n, 

„>^ + ^|l + l ,l + ^^]. (24, 

In this way, convergence of Gff(E) can be proven for arbitrary model parameters. (An alternative would be to apply 
the corollary used by Schweber [31]). 

I shall now demonstrate that this proof of convergence is not sufficient to show that lini/v->oc Gq N \e) is related 
to the spectrum of the full Rabi model for fixed parity. The reason is the intrinsic ambiguity in the definition of the 
finite-dimensional approximants to Hr. If Hr were compact, it could be defined as the limit limjv->oo H^ N > in the 
norm topology [29] and a(H^) would converge towards o{Hr) in the sense defined in (13) because the sets {E^} 
and {E n } would satisfy the even stronger condition 

lim sup \EW-E n \ = 0. (25) 

W^oo„ eINo 

But Hr is unbounded in H and therefore not the limit of any sequence of finite-dimensional operators. Even if the 
sequence of truncated state spaces {"H^jneiN could be determined in a unique way (which is not the case), the 
approximating operators need not to be the H^ — P^HrP^ n \ used in method B. The projection pW 

affects only the states \N) and \N + 1) because it is equivalent to setting (N\Hr\N + 1) = \g\VN + 1 to zero. But 
as these matrix elements grow as \/N for N — > oo, it is not clear whether the spectra a(H^) converge to <j(Hr), 
even if they do converge to a limit. One could think of other modifications of the matrix elements of Hr between 
states with high energy and claim that these are also admissible truncation prescriptions. A similar ambiguity is 
well-known from quantum field theory: The regularization scheme needed to make the perturbative expressions finite 
by truncating the ultra-violet modes is not uniquely definable. One has to prove that the resulting (renormalized) 
expressions are independent from the choice of cut-off. In our case, the task is to show that the spectra of a sequence 
of finite-dimensional Hamiltonians {H^} converge to o{Hr), and to this end one would have to prove that this limit 
does not depend on the modification of Hr in the high-energy region which is used to define the approximants . 
The formal condition reads as follows: 

If the elements of two sequences of operators {J?W} and {H^} 

are defined on the spaces H (N) for N = 0,1,... 
with H (n) C n w C . . .L 2 (K) and = = H R if projected onto H^-^, the limit of the spectra of both 

sequences as defined in (13) must coincide: limAr^oo a{H^) = liniTv^oo cr(H^). 

This is obviously necessary (albeit not sufficient) to justify the claim that a(H^) converges to the unique <j(Hr), 
if Hr is not compact and the H W do not converge to it in the operator sense, but it is violated in the present case: 
One may define a cut-off prescription for H( N ) which affects only the states | N — 1), | N) and \N + 1) but produces a 

(fictitious) pole in Gq N \e) at an arbitrary value E a . Consider the expression (20) for Go, which reads 

G Q (E) = l\b (E) - ai \h(E) - a 2 \b 2 (E) . . . a N ^\b N ^ - a N G N (E) (26) 

as follows from (19). The standard cut-off using the projection pW sets G N (E) = l/b N (E). But G N {E) can also 
be determined from below by a recurrence relation equivalent to (19), 

^^-^-^TjE-y (27) 

a,j+i a,j + iLTj{Ji,) 

If one sets G (E ) = oo, corresponding to a pole of G {E) at E a , the recurrence (27) with the initial condition 
Gi(^o) = bQ{E )/ai yields a G N (E ) which will tend to the value -uj/g 2 for N -> oo. The matrix clement (N\H^\N) 
of the truncated Hamiltonian is defined by 

(N\HW\N) = H%% =Eo- ^-y. (28) 
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For all 7^ (N, N), H^ — H^\ But the ensuing Gq N \e) will have a pole at Eq, although only a single matrix 
element was changed and therefore _p( Ar_1 )H"( w )p( Ar_1 ) = p( N - x ) HrP^ n ~ x \ The same result follows if one sets 
^nn = Eo — N/Gn(Eo) and modifies at the same time H]$_i n ~ 9^ 1 which does not yield a low energy state a 

priori if GN+oj/g 2 tends slower to zero than 1/N. It is therefore possible to define a sequence of Hamiltonians {H^} 
with eigenvalues which are not close to elements of the spectrum of the untruncated Rabi model, but differ from the 
H( N ~) in at most three matrix elements and coincide with Hr on 'H^ N ~ 1 \ Neither {H^} nor {ifW} converge to Hr 
in norm topology because Hr is unbounded. This entails that no independent criterion exists to select the "correct" 
sequence {H^} apart from convergence of the spectrum in the sense of (13) to some set {e„}„ e iN . However, the 
spectra of the sequence {H^} converge as well, if H^^, H^l f _ 1 are appropriate functions of N. The problem 
cannot be solved by consideration of resolvents instead of the Hamiltonians. Whereas Hr has compact resolvent, 
the resolvents of the finite-range operators ifW and H^ are bounded but not compact and do not converge to 
(z — Hr)^ 1 in norm topology. Even if one could prove that only {ffW} and not {H^} converges to Hr in strong 
resolvent sense, it would be insufficient to determine ct(Hr) unambiguously [30]. 

The standard cut-off prescription leads indeed to the correct spectrum of the full model as the comparison with 
the analytical solution [3] shows. But this prescription is only valid because the finite approximants Hn obtained 
by method A are the projections PnHrPn and thus correspond to the standard cut-off of the model in the finite- 
dimensional state spaces 'Hat as shown above. These spaces are not generated by the first N + 1 Fock states as the 
%W but are sufficiently close to them for large ./V because the coherent factor e~ gz ^ appearing in the ip(z; N) has a 
convergent expansion in Fock states. The isomorphism I between B and i 2 (R) maps the function z" j\fn\ onto |n); 
it follows that for any vector \iPn) G I(Hn) C L 2 (R) and its projection \ip' N ) onto %^ N \ we have 

lim (Vjv - V n \>I>n ~ i>' N ) = 0- (29) 

Therefore P N H R P N - P^HrP^ tends to zero in the norm topology for N — > 00. We may conclude that method 
B yields the spectrum of the Rabi model in the limit N — > 00 precisely because the finite-dimensional approximants 
H N of method A are the correct ones. Method B cannot be justified independently from method A [32]. 

Apart from the implicit dependence on method A, method B shares with it the same problems: The expressions 
derived from continued fractions do not allow for qualitative analysis because the singularities of neither F^E) 
nor 1/Gq(E) are known, which would be necessary to infer the distribution of the eigenvalues of Hr. Moreover, 
the position and nature of spectral degeneracies cannot be obtained from those functions, as exemplified in [14], 
where the following is said on level crossings: "The individual matrix elements g, h avoid level crossing due to parity 
conservation, since eigenstates of consecutive eigenvalues have different parity." However, g and h correspond to the 
subspaces with even and odd parity, respectively; the parity does not change within each subspace. The idea that 
alternating parity forbids level crossing between consecutive eigenvalues is erroneous, because the converse is true [3]. 
If applied to the full Hilbcrt space, Zieglers statement is wrong as well: Energetically neighboring states do not have 
always different parity. Finally, the fact that the degeneracies (of states with different parity) occur only if E + g 2 /u> 
is an integer multiple of uj can be deduced neither with method A nor method B. 



IV. CONCLUSIONS 



We conclude that both methods A and B yield correct approximations for the spectrum of the quantum Rabi model 
and are equivalent to numerical exact diagonalization in finite-dimensional Hilbert spaces of sufficient large dimension. 
Diagonalization on a truncated state space has been used for many years as a tool to obtain numerically exact results 
which were subsequently compared to a variety of analytical approximations (see e.g. [33]) without questioning 
its correctness, although the Hamiltonian Hr is not the limit of a sequence of finite-dimensional operators. The 
soundness of the approach was taken for granted because of numerical convergence of the results, which is equivalent 
with convergence of the corresponding continued fraction. The pointwise convergence of the latter has been established 
by Schweber [5] for method A and by Durst et al. [26] for method B. In the case of method A, this convergence is 
equivalent to the analyticity condition for the untruncated model and therefore sufficient to prove numerical identity 
with the exact spectrum for large enough N. In contrast to method A, the convergence of the continued fraction says 
nothing about the spectrum of the full model within the framework of method B. The sequence of approximating 
Hamiltonians H^ N ^ is not uniquely determined and spectral convergence alone is not sufficient to prove a connection 
with Hr. But because the cut-off scheme employed by method B is equivalent (for large N) to the sequence of 
finite-dimensional approximants Hm from method A, it yields the (numerically) exact spectrum in the limit N — > 00 
as well. It owes this to the validity of method A: Contrary to the claim made in [14], method B cannot be considered 
as an independent way to obtain the correct spectrum of Hr. 
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In this way the use of continued fractions and (equivalently) numerical exact diagonalization can be rigorously 
justified in the case of the Rabi model. The key to the proof is Bargmanns analyticity condition for functions 
as elements of an infinite-dimensional Hilbert space. Whether this equivalence with the analytical solution can be 
extended to more complicated systems or breaks down in some cases [34] should be the subject of future study. 
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